### PPP
### Mediation, VaR, 

rm(list = ls())

### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

### libraries and packages
pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", "grid",
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", "gridExtra",
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)

load("ppp_cleaned")

p2 <- as.data.frame(cbind(p2$l_sat, p2$c_sat, p2$l_watch, p2$c_watch, p2$econ, p2$anti_west, p2$day_pollute, 
                          p2$life_sat,
                          p2$day, p2$parade, p2$aqi, p2$aqi.lagged, p2$aqi.lagged2, p2$national, p2$hukou, 
                    p2$age, p2$income, p2$gender, p2$ccp, p2$soe, p2$insider, p2$married, p2$educ))
colnames(p2) <- c("l_sat", "c_sat", "l_watch", "c_watch", "econ", "anti_west", "day_pollute", 
                  "life_sat","day", "parade", "aqi", "aqi.lagged", "aqi.lagged2", "national", "hukou", 
                        "age", "income", "gender", "ccp", "soe", "insider", "married", "educ")
p2 <- na.omit(p2)

###########
setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Figures/Mediation")
set.seed(123)

med.fit <- lm(aqi ~ parade + aqi.lagged + aqi.lagged2 + national + gender + ccp + hukou + age + income + gender + ccp + soe + insider + married + educ, data = p2)
out.fit <- lm(l_sat ~ parade + aqi.lagged + aqi.lagged2 + national + aqi + gender + ccp + hukou + age + income + gender + ccp + educ + married + soe + insider, data = p2)
med.out <- mediate(med.fit, out.fit, treat = "parade", mediator = "aqi",  dropobs = TRUE, sims = 1000)

pdf(file = "l_sat_mediation.pdf", height = 6, width = 10)
plot(med.out, main = "Through reducing pollution, the parade \n enhances citizens' satisfaction with the local government")
dev.off()
# ACME
ACME_lsat <- med.out$d0
ACME_lsat_lc <- med.out$d0.ci[1]
ACME_lsat_uc <- med.out$d0.ci[2]
ADE_lsat <- med.out$z0
ADE_lsat_lc <- med.out$z0.ci[1]
ADE_lsat_uc <- med.out$z0.ci[2]
TE_lsat <- med.out$tau.coef
TE_lsat_lc <- med.out$tau.ci[1]
TE_lsat_uc <- med.out$tau.ci[2]

# c_sat
med.fit <- lm(aqi ~ parade + aqi.lagged + aqi.lagged2 + national + gender + ccp + hukou + age + income + gender + ccp + soe + insider + married + educ, data = p2)
out.fit <- lm(c_sat ~ parade + aqi.lagged + aqi.lagged2 + national + aqi + gender + ccp + hukou + age + income + gender + ccp + educ + married + soe + insider, data = p2)
med.out <- mediate(med.fit, out.fit, treat = "parade", mediator = "aqi",
                     dropobs = TRUE, sims = 1000)
pdf(file = "c_sat_mediation.pdf", height = 6, width = 10)
plot(med.out, main = "Through reducing pollution, the parade \n enhances citizens' satisfaction with the central government")
dev.off()
summary(med.out)
ACME_csat <- med.out$d0
ACME_csat_lc <- med.out$d0.ci[1]
ACME_csat_uc <- med.out$d0.ci[2]
ADE_csat <- med.out$z0
ADE_csat_lc <- med.out$z0.ci[1]
ADE_csat_uc <- med.out$z0.ci[2]
TE_csat <- med.out$tau.coef
TE_csat_lc <- med.out$tau.ci[1]
TE_csat_uc <- med.out$tau.ci[2]

# l_watchdog
med.fit <- lm(aqi ~ parade + national + aqi.lagged + aqi.lagged2 + gender + ccp + hukou + age + income + gender + ccp + soe + insider + married + educ, data = p2)
out.fit <- lm(l_watch ~ parade + aqi.lagged + aqi.lagged2 + aqi + gender + ccp + national + hukou + age + income + gender + ccp + educ + married + soe + insider, data = p2)
med.out <- mediate(med.fit, out.fit, treat = "parade", mediator = "aqi",
                     dropobs = TRUE, sims = 1000)
pdf(file = "l_watchdog_mediation.pdf", height = 6, width = 10)
plot(med.out, main = "Through reducing pollution, the parade \n dampens citizens' need 
     for watchdogs at the local level")
dev.off()
summary(med.out)

ACME_lwatch <- med.out$d0
ACME_lwatch_lc <- med.out$d0.ci[1]
ACME_lwatch_uc <- med.out$d0.ci[2]
ADE_lwatch <- med.out$z0
ADE_lwatch_lc <- med.out$z0.ci[1]
ADE_lwatch_uc <- med.out$z0.ci[2]
TE_lwatch <- med.out$tau.coef
TE_lwatch_lc <- med.out$tau.ci[1]
TE_lwatch_uc <- med.out$tau.ci[2]



# c_watchdog
med.fit <- lm(aqi ~ parade + aqi.lagged + aqi.lagged2 + national + gender + ccp + hukou + age + income + gender + ccp + soe + insider + married + educ, data = p2)
out.fit <- lm(c_watch ~ parade + aqi.lagged + aqi.lagged2 + national + aqi + gender + ccp + hukou + age + income + gender + ccp + educ + married + soe + insider, data = p2)
med.out <- mediate(med.fit, out.fit, treat = "parade", mediator = "aqi",
                     dropobs = TRUE, sims = 1000)
pdf(file = "c_watchdog_mediation.pdf", height = 6, width = 10)
plot(med.out, main = "Through reducing pollution, the parade \n dampens citizens' need 
     for watchdogs at the central level")
dev.off()
summary(med.out)

ACME_cwatch <- med.out$d0
ACME_cwatch_lc <- med.out$d0.ci[1]
ACME_cwatch_uc <- med.out$d0.ci[2]
ADE_cwatch <- med.out$z0
ADE_cwatch_lc <- med.out$z0.ci[1]
ADE_cwatch_uc <- med.out$z0.ci[2]
TE_cwatch <- med.out$tau.coef
TE_cwatch_lc <- med.out$tau.ci[1]
TE_cwatch_uc <- med.out$tau.ci[2]

# Anti-west
med.fit <- lm(aqi ~ parade + aqi.lagged + aqi.lagged2 + national + gender + ccp + hukou + age + income + gender + ccp + soe + insider + married + educ, data = p2)
out.fit <- lm(anti_west ~ parade + aqi.lagged + aqi.lagged2 + national + aqi + gender + ccp + hukou + age + income + gender + ccp + educ + married + soe + insider, data = p2)
med.out <- mediate(med.fit, out.fit, treat = "parade", mediator = "aqi",
                     dropobs = TRUE, sims = 1000)
pdf(file = "anti_west_mediation.pdf", height = 6, width = 10)
plot(med.out, main = "Does the parade fuel anti-west sentiment through reducing pollution?")
dev.off()
summary(med.out)

ACME_antiwest <- med.out$d0
ACME_antiwest_lc <- med.out$d0.ci[1]
ACME_antiwest_uc <- med.out$d0.ci[2]
ADE_antiwest <- med.out$z0
ADE_antiwest_lc <- med.out$z0.ci[1]
ADE_antiwest_uc <- med.out$z0.ci[2]
TE_antiwest <- med.out$tau.coef
TE_antiwest_lc <- med.out$tau.ci[1]
TE_antiwest_uc <- med.out$tau.ci[2]


# subjective well-being  
med.fit <- lm(aqi ~ parade + aqi.lagged + aqi.lagged2 + national + gender + ccp + hukou + age + income + gender + ccp + soe + insider + married + educ, data = p2)
out.fit <- lm(life_sat ~ parade + aqi.lagged + aqi.lagged2 + national + aqi + gender + ccp + hukou + age + income + gender + ccp + educ + married + soe + insider, data = p2)
med.out <- mediate(med.fit, out.fit, treat = "parade", mediator = "aqi",
                     dropobs = TRUE, sims = 1000)
pdf(file = "life_sat_mediation.pdf", height = 6, width = 10)
plot(med.out, main = "Parade enhances life satisfaction through reducing actual pollution")
dev.off()
summary(med.out)

ACME_life <- med.out$d0
ACME_life_lc <- med.out$d0.ci[1]
ACME_life_uc <- med.out$d0.ci[2]
ADE_life <- med.out$z0
ADE_life_lc <- med.out$z0.ci[1]
ADE_life_uc <- med.out$z0.ci[2]
TE_life <- med.out$tau.coef
TE_life_lc <- med.out$tau.ci[1]
TE_life_uc <- med.out$tau.ci[2]

###### Making seven figures with ggplot2 ######
## lsat ##
coef.vect <- c(ACME_lsat, ADE_lsat, TE_lsat)
lower <- c(ACME_lsat_lc, ADE_lsat_lc,TE_lsat_lc)
upper <- c(ACME_lsat_uc, ADE_lsat_uc, TE_lsat_uc)

longnames <- factor(c("Mediation effect of \n Pollution reduction", "Direct effect of the parade", 
                      "Total effect of the parade"),
                    levels=c("Total effect of the parade", "Direct effect of the parade", 
                             "Mediation effect of \n Pollution reduction"))

frame = data.frame(longnames, coef.vect, upper, lower)
c1 <- ggplot(frame, aes(longnames, coef.vect))
d1 <- c1 + geom_point(size=3) + geom_errorbar(aes(ymin=lower,
                                           ymax=upper, width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("") + ylab("Coefficients") +
  ggtitle("Local Gov. Satisfaction") +
  theme(plot.title = element_text(vjust = 1.3, hjust=.5))

## csat ##
coef.vect <- c(ACME_csat, ADE_csat, TE_csat)
lower <- c(ACME_csat_lc, ADE_csat_lc,TE_csat_lc)
upper <- c(ACME_csat_uc, ADE_csat_uc, TE_csat_uc)

longnames <- factor(c("Mediation effect of \n Pollution reduction", "Direct effect of the parade", 
                      "Total effect of the parade"),
                    levels=c("Total effect of the parade", "Direct effect of the parade", 
                             "Mediation effect of \n Pollution reduction"))

frame = data.frame(longnames, coef.vect, upper, lower)
c2 <- ggplot(frame, aes(longnames, coef.vect))
d2 <- c2 + geom_point(size=3) + geom_errorbar(aes(ymin=lower,
                                           ymax=upper, width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("") + ylab("Coefficients") +
  ggtitle("Central Gov. Satisfaction") +
  theme(plot.title = element_text(vjust = 1.3, hjust=.5))

## lwatch ##
coef.vect <- c(ACME_lwatch, ADE_lwatch, TE_lwatch)
lower <- c(ACME_lwatch_lc, ADE_lwatch_lc,TE_lwatch_lc)
upper <- c(ACME_lwatch_uc, ADE_lwatch_uc, TE_lwatch_uc)

longnames <- factor(c("Mediation effect of \n Pollution reduction", "Direct effect of the parade", 
                      "Total effect of the parade"),
                    levels=c("Total effect of the parade", "Direct effect of the parade", 
                             "Mediation effect of \n Pollution reduction"))

frame = data.frame(longnames, coef.vect, upper, lower)
c3 <- ggplot(frame, aes(longnames, coef.vect))
d3 <- c3 + geom_point(size=3) + geom_errorbar(aes(ymin=lower,
                                           ymax=upper, width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("") + ylab("Coefficients") +
  ggtitle("Local Gov. Oversight") +
  theme(plot.title = element_text(vjust = 1.3, hjust=.5), title = element_text(size = 10))

## cwatch ##
coef.vect <- c(ACME_cwatch, ADE_cwatch, TE_cwatch)
lower <- c(ACME_cwatch_lc, ADE_cwatch_lc,TE_cwatch_lc)
upper <- c(ACME_cwatch_uc, ADE_cwatch_uc, TE_cwatch_uc)

longnames <- factor(c("Mediation effect of \n Pollution reduction", "Direct effect of the parade", 
                      "Total effect of the parade"),
                    levels=c("Total effect of the parade", "Direct effect of the parade", 
                             "Mediation effect of \n Pollution reduction"))

frame = data.frame(longnames, coef.vect, upper, lower)
c4 <- ggplot(frame, aes(longnames, coef.vect))
d4 <- c4 + geom_point(size=3) + geom_errorbar(aes(ymin=lower,
                                           ymax=upper, width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("") + ylab("Coefficients") +
  ggtitle("Central Gov. Oversight") +
  theme(plot.title = element_text(vjust = 1.3, hjust=.5))


## anti-west ##
coef.vect <- c(ACME_antiwest, ADE_antiwest, TE_antiwest)
lower <- c(ACME_antiwest_lc, ADE_antiwest_lc,TE_antiwest_lc)
upper <- c(ACME_antiwest_uc, ADE_antiwest_uc, TE_antiwest_uc)

longnames <- factor(c("Mediation effect of \n Pollution reduction", "Direct effect of the parade", 
                      "Total effect of the parade"),
                    levels=c("Total effect of the parade", "Direct effect of the parade", 
                             "Mediation effect of \n Pollution reduction"))

frame = data.frame(longnames, coef.vect, upper, lower)
c6 <- ggplot(frame, aes(longnames, coef.vect))
d6 <- c6 + geom_point(size=3) + geom_errorbar(aes(ymin=lower,
                                           ymax=upper, width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("") + ylab("Coefficients") +
  ggtitle("Anti-West") +
  theme(plot.title = element_text(vjust = 1.3, hjust=.5),title = element_text(size = 10.5))

## LIFE ##
coef.vect <- c(ACME_life, ADE_life, TE_life)
lower <- c(ACME_life_lc, ADE_life_lc,TE_life_lc)
upper <- c(ACME_life_uc, ADE_life_uc, TE_life_uc)

longnames <- factor(c("Mediation effect of \n Pollution reduction", "Direct effect of the parade", 
                      "Total effect of the parade"),
                    levels=c("Total effect of the parade", "Direct effect of the parade", 
                             "Mediation effect of \n Pollution reduction"))

frame = data.frame(longnames, coef.vect, upper, lower)
c7 <- ggplot(frame, aes(longnames, coef.vect))
d7 <- c7 + geom_point(size=3) + geom_errorbar(aes(ymin=lower,
                                           ymax=upper, width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("") + ylab("Coefficients") +
  ggtitle("Life satisfaction") +
  theme(plot.title = element_text(vjust = 1.3))


setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Figures/Mediation")

pdf(file = "mediation_in_one_landscape.pdf", height = 4, width = 12)
grid.arrange(arrangeGrob(d1,
                         d3,
                         d2,
                         d4,
                         d6,
                         ncol=3), heights=c(13, 1))
dev.off()

## TABLE A10 
pdf(file = "mediation_in_one_portrait.pdf", height = 6, width = 8)
grid.arrange(arrangeGrob(d1,
                         d2,
                         d3,
                         d4,
                         d6,
                         ncol=2), heights=c(13, 1))
dev.off()